load libraries

library(tidyverse)
library(psych)
library(knitr)
library(lme4)
library(lmerTest)
devtools::install_github('jflournoy/qualtrics')

set color palettes

palette = wesanderson::wes_palette("Zissou1", n = 4, type = "continuous")
palette_cond = wesanderson::wes_palette("Darjeeling1", 3, "continuous")
palette_grey = "#54565B"

load data

# define variables and paths
sub_dir = "~/Dropbox (PfeiBer Lab)/FreshmanProject/Tasks/ROC-C/output/"
sub_pattern = "HL[0-9]{3}"
subjects = list.files(sub_dir, pattern = sub_pattern)
runs = c("run1", "run2", "run3")

# initialize data frame
data = data.frame()

# load data
for (sub in subjects) {
  for (run in runs) {
    file = paste0(sub_dir, sub, '/', sub, '_', run,'.csv')
    tmp = tryCatch(read.csv(file, stringsAsFactors = FALSE) %>%
                   mutate(subjectID = sub,
                          run = run,
                          respCue = as.integer(as.character(respCue)),
                          respRating = as.integer(as.character(respRating)),
                          respEffort = as.integer(as.character(respEffort))), error = function(e) NULL)
      data = bind_rows(data, tmp)
      rm(tmp)
  }
}

# load condition info
condition = read.csv("Himmel Info - Himmel III Participant Info.csv", stringsAsFactors = FALSE) %>%
  rename("subjectID" = Subject.ID,
         "condition" = Condition,
         "experimenter" = Experimenter) %>%
  mutate(condition = ifelse(condition == "A", "control",
                     ifelse(condition == "B", "food",
                     ifelse(condition == "C", "autonomy", NA)))) %>%
  select(subjectID, condition, experimenter)

get experimental manipulation survey data

# define variables
cred_file_location = '~/credentials.yaml.DEFAULT'
sid_column_name = '(SID)'
survey_name_filter = 'Autonomy Manipulation'
sid_pattern = 'HL[0-9]{3}'
exclude_sid = 'HL999' # subject IDs to exclude
identifiableData = c('IPAddress') # exclude when printing duplicates

# load credential file
credentials = scorequaltrics::creds_from_file(cred_file_location)

# filter
surveysAvail = scorequaltrics::get_surveys(credentials)
surveysFiltered = filter(surveysAvail, grepl(survey_name_filter, SurveyName))

# get survey
surveys = scorequaltrics::get_survey_responses(credentials, 
                                               surveyid = surveysFiltered$SurveyID[[1]])

# tidy surveys
surveys1 = surveys %>%
  slice(-1:-48) %>% # filter out pilot observations in the first 48 rows 
  mutate(subjectID = ifelse(grepl("^[0-9]{3}", SID), paste0("HL", SID),
                     ifelse(grepl("HL[0-9]{3}", SID), SID, NA))) %>%
  filter(!subjectID == "HL999")

# check number of observations
surveys1 %>%
  group_by(subjectID) %>%
  summarize(n = n()) %>%
  arrange(desc(n))
# filter our duplicate subject and recode condition B choices
manipulation = surveys1 %>%
  filter(!qid == "R_3HC8QsQwOBNYF9z") %>% # filter out first observation of HL009
  select(subjectID, COND, B_1) %>%
  rename("condition" = "COND") %>%
  mutate(condition = ifelse(condition == "1", "food",
                     ifelse(condition == "2", "autonomy", NA)),
         taskChoice = ifelse(B_1 == "1", "regulate",
                   ifelse(B_1 == "2", "look", NA))) %>%
  select(-B_1)

get post-task survey data

# define variables
survey_name_filter = 'Himmel Survey 2018'

# filter
surveysAvail = scorequaltrics::get_surveys(credentials)
surveysFiltered = filter(surveysAvail, grepl(survey_name_filter, SurveyName))

# get survey
survey.post = scorequaltrics::get_survey_responses(credentials, 
                                               surveyid = surveysFiltered$SurveyID[[1]])

# tidy surveys
survey.post1 = survey.post %>%
  mutate(subjectID = ifelse(grepl("^[0-9]{3}", SID), paste0("HL", SID),
                     ifelse(grepl("HL[0-9]{3}", SID), SID, NA))) %>%
  filter(!subjectID == "HL999")

# check number of observations
survey.post1 %>%
  group_by(subjectID) %>%
  summarize(n = n()) %>%
  arrange(desc(n))
# filter our duplicate subject and recode condition B choices
check = survey.post1 %>%
  filter(!qid == "R_3OqwT3EzxpDlfOI") %>% # filter out first observation of HL009
  select(subjectID, starts_with("CHECK"), starts_with("VALUES")) %>%
  gather(set, value, starts_with("CHECK")) %>%
  extract(set, c("set", "question"), "(CHECK_[1-2]{1})_([1-4]{1})") %>%
  mutate(choice = ifelse(set == "CHECK_1", "yes",
                  ifelse(set == "CHECK_2", "no", NA)),
         question = ifelse(question == "1", "liking",
                    ifelse(question == "2", "engagement",
                    ifelse(question == "3", "motivation",
                    ifelse(question == "4", "difficulty", NA)))),
         value = as.numeric(value)) %>%
  select(-set, -starts_with("VALUES"))

values = survey.post1 %>%
  filter(!qid == "R_3OqwT3EzxpDlfOI") %>% # filter out first observation of HL009
  select(subjectID, starts_with("VALUES")) %>%
  rename("Being healthy" = VALUES_1,
         "Being physically fit" = VALUES_2,
         "Eating healthfully" = VALUES_3,
         "Making decisions for myself" = VALUES_4,
         "Having a choice" = VALUES_5,
         "Being independent" = VALUES_6,
         "Enjoying what I eat" = VALUES_7,
         "Controlling my desire for unhealthy foods" = VALUES_8,
         "Not having to make decisions" = VALUES_9,
         "Getting help from others to make decisions" = VALUES_10,
         "Minimizing the number of decisions I make" = VALUES_11,
         "Following instructions" = VALUES_12) %>%
  gather(question, value, -subjectID)

IMI = survey.post1 %>%
  filter(!qid == "R_3OqwT3EzxpDlfOI") %>% # filter out first observation of HL009
  select(subjectID, starts_with("IMI")) %>%
  gather(item, value, -subjectID) %>%
  mutate(factor = ifelse(item %in% sprintf("IMI_%d", 1:7), "interest_enjoyment",
                  ifelse(item %in% sprintf("IMI_%d", 8:13), "competence",
                  ifelse(item %in% sprintf("IMI_%d", 14:18), "effort_importance",
                  ifelse(item %in% sprintf("IMI_%d", 19:23), "pressure_tension",
                  ifelse(item %in% sprintf("IMI_%d", 24:30), "perceived_choice",
                  ifelse(item %in% sprintf("IMI_%d", 31:37), "value_usefulness", NA)))))),
         value = ifelse(item %in% sprintf("IMI_%d", c(3,4,13,15,18,19,21,25,26,27,28,30)) & value == 1, 7,
                 ifelse(item %in% sprintf("IMI_%d", c(3,4,13,15,18,19,21,25,26,27,28,30)) & value == 2, 6,
                 ifelse(item %in% sprintf("IMI_%d", c(3,4,13,15,18,19,21,25,26,27,28,30)) & value == 3, 5,
                 ifelse(item %in% sprintf("IMI_%d", c(3,4,13,15,18,19,21,25,26,27,28,30)) & value == 5, 3,
                 ifelse(item %in% sprintf("IMI_%d", c(3,4,13,15,18,19,21,25,26,27,28,30)) & value == 6, 2,
                 ifelse(item %in% sprintf("IMI_%d", c(3,4,13,15,18,19,21,25,26,27,28,30)) & value == 7, 1, value))))))) %>%
  group_by(subjectID, factor) %>%
  summarize(mean = mean(as.numeric(value), na.rm = TRUE))

tidy data

data1 = data %>%
  filter(!subjectID %in% c("HL999")) %>%
  mutate(rtCue = ifelse(rtCue == "NaN", NA, rtCue),
         rtRating = ifelse(rtRating == "NaN", NA, rtRating),
         rtEffort = ifelse(rtEffort == "NaN", NA, rtEffort),
         action = ifelse(respCue == 1, "look",
                  ifelse(respCue == 2, "regulate", NA)),
         action = ifelse(cond == "LOOK" & is.na(respCue), "look",
                  ifelse(cond == "REGULATE" & is.na(respCue), "regulate", action)),
         action = as.factor(action),
         choice = ifelse(cond %in% c("LOOK", "REGULATE"), "no",
                  ifelse(cond == "CHOOSE", "yes", NA)),
         choice = as.factor(choice),
         ratingDiff = respRating - likingRating) %>%
  extract(col = foodPic, into = "category", regex = "(.*)[0-9]{2}", remove = FALSE) %>%
  group_by(subjectID) %>%
  mutate(trial = row_number()) %>%
  left_join(condition, by = "subjectID") %>%
  left_join(manipulation, by = c("subjectID", "condition")) %>%
  left_join(check, by = c("subjectID", "choice")) %>%
  select(subjectID, condition, taskChoice, run, action, choice, cond, respCue, trial, everything()) %>%
  group_by(subjectID) %>%
  mutate(meanLiking = mean(likingRating, na.rm = TRUE))
## Warning: Column `choice` joining factor and character vector, coercing into
## character vector

engagement by task and experimental condition

data1 %>%
  filter(!is.na(question)) %>%
  group_by(subjectID, condition, choice) %>%
  ggplot(aes(condition, value, color = choice)) +
    stat_summary(aes(group = choice), fun.y = mean, geom = "line") +
    stat_summary(fun.data = "mean_cl_boot") +
    facet_grid(~question) + 
    scale_color_manual(values = c(palette[1], palette[3])) + 
    theme_minimal()

intrinsic motivation

IMI %>%
  left_join(., select(data1, subjectID, condition), by = "subjectID") %>%
  ggplot(aes(condition, mean, color = factor)) +
    stat_summary(aes(group = factor), fun.y = mean, geom = "line") +
    stat_summary(fun.data = "mean_cl_boot", size = .5) +
    scale_color_manual(values = wesanderson::wes_palette("Zissou1", 7, "continuous")) + 
    theme_minimal()

IMI %>%
  left_join(., select(data1, subjectID, condition), by = "subjectID") %>%
  ggplot(aes(factor, mean, color = condition)) +
    stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
    stat_summary(fun.data = "mean_cl_boot", size = .5) +
    scale_color_manual(values = wesanderson::wes_palette("Darjeeling1", 3, "continuous")) + 
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

values %>%
  left_join(., select(data1, subjectID, condition), by = "subjectID") %>%
  ggplot(aes(condition, value, color = question)) +
    stat_summary(aes(group = question), fun.y = mean, geom = "line") +
    stat_summary(aes(group = question), fun.data = "mean_cl_boot", size = .5) +
    scale_color_manual(values = wesanderson::wes_palette("Zissou1", 12, "continuous")) + 
    theme_minimal()

values %>%
  left_join(., select(data1, subjectID, condition), by = "subjectID") %>%
  ggplot(aes(question, value, color = condition)) +
    stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
    stat_summary(aes(group = condition), fun.data = "mean_cl_boot", size = .5) +
    scale_color_manual(values = wesanderson::wes_palette("Darjeeling1", 3, "continuous")) + 
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

missing choice data

missing.plot = data1 %>%
  filter(choice == "yes") %>%
  group_by(subjectID, action) %>%
  summarize(n = n()) %>%
  spread(action, n) %>%
  rename("missing" = `<NA>`) %>%
  mutate(missing = ifelse(is.na(missing), 0, missing),
         missingSets = missing / 3,
         percentMissing = (missing / (look + regulate + missing)) * 100)

ggplot(missing.plot, aes(missingSets)) +
  geom_histogram(bins = 4) +
  labs(x = "number of sets missing") + 
  theme_minimal()

missing.plot %>%
  ungroup() %>%
  summarize(n = n(),
            mean = mean(percentMissing, na.rm = TRUE),
            sd = sd(percentMissing, na.rm = TRUE),
            min = min(percentMissing, na.rm = TRUE),
            max = max(percentMissing, na.rm = TRUE)) %>%
      kable(digits = 2, format = "pandoc", caption = "average percent of choice trials missing")
average percent of choice trials missing
n mean sd min max
105 3.02 4.56 0 16.67

percent regulate

all conditions

There are a few outliers, but most people chose to regulate around 50% of the time.

percent.plot = data1 %>%
  filter(choice == "yes") %>%
  group_by(subjectID, action, condition, taskChoice) %>%
  summarize(n = n()) %>%
  spread(action, n) %>%
  mutate(percent.regulate = (regulate / (look + regulate)) * 100) %>%
  arrange(percent.regulate)

ggplot(percent.plot, aes("", percent.regulate)) +
  geom_boxplot() +
  facet_grid(~condition) + 
  labs(x = "", y = "percent chosen to regulate") + 
  theme_minimal()

percent.plot %>%
  group_by(condition) %>%
  summarize(n = n(),
            mean = mean(percent.regulate, na.rm = TRUE),
            sd = sd(percent.regulate, na.rm = TRUE),
            min = min(percent.regulate, na.rm = TRUE),
            max = max(percent.regulate, na.rm = TRUE)) %>%
      kable(digits = 2, format = "pandoc", caption = "group average by condition")
group average by condition
condition n mean sd min max
autonomy 32 50.34 9.07 29.41 73.33
control 35 49.75 6.95 27.78 66.67
food 38 47.36 8.69 25.00 66.67
percent.plot

by profile in food manipulation condition

ggplot(filter(percent.plot, condition == "food"), aes("", percent.regulate)) +
  geom_boxplot() +
  facet_grid(~taskChoice) + 
  labs(x = "", y = "percent chosen to regulate", title = "percent regulation choices in food manipulation condition only") + 
  theme_minimal()

percent.plot %>%
  group_by(condition, taskChoice) %>%
  summarize(n = n(),
            mean = mean(percent.regulate, na.rm = TRUE),
            sd = sd(percent.regulate, na.rm = TRUE),
            min = min(percent.regulate, na.rm = TRUE),
            max = max(percent.regulate, na.rm = TRUE)) %>%
      kable(digits = 2, format = "pandoc", caption = "group average by condition")
group average by condition
condition taskChoice n mean sd min max
autonomy NA 32 50.34 9.07 29.41 73.33
control NA 35 49.75 6.95 27.78 66.67
food look 22 48.39 8.64 25.00 66.67
food regulate 16 45.95 8.84 33.33 58.82

liking descriptives

mean liking ratings by subject

data1 %>%
  group_by(subjectID) %>%
  summarize(mean = mean(meanLiking, na.rm = TRUE)) %>%
  ggplot(aes(mean)) +
    geom_histogram(stat = "bin", bins = 10) +
    theme_minimal()

data1 %>%
  group_by(subjectID) %>%
  summarize(mean = mean(meanLiking, na.rm = TRUE)) %>%
  mutate(lt2.5 = ifelse(mean < 2.5, 1, 0),
         lt3 = ifelse(mean < 3, 1, 0)) %>%
  summarize(`less than 2.5` = sum(lt2.5),
            `less than 3` = sum(lt3))

craving ratings

by subject and task condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition) %>%
  summarize(meanRating = mean(respRating, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanRating)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = c(palette_cond)) +
  labs(y = "mean craving rating") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanRating)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = c(palette_cond, palette[1],palette[3])) +
  labs(y = "mean craving rating") + 
  theme_minimal()

data1 %>%
  group_by(action, choice) %>%
  summarize(meanRating = mean(respRating, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean craving ratings by condition")
mean craving ratings by condition
action choice meanRating n
look no 2.96 7452
look yes 2.99 11061
regulate no 2.35 7428
regulate yes 2.31 10608
NA yes 2.80 651

by task and experimental condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition, taskChoice) %>%
  summarize(meanRating = mean(respRating, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean craving rating") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean craving rating") + 
  theme_minimal()

data1 %>%
  group_by(action, choice, condition) %>%
  summarize(meanRating = mean(respRating, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean craving ratings by condition")
mean craving ratings by condition
action choice condition meanRating n
look no autonomy 3.04 2262
look no control 2.97 2454
look no food 2.91 2736
look yes autonomy 3.05 3243
look yes control 2.99 3630
look yes food 2.95 4188
regulate no autonomy 2.40 2238
regulate no control 2.26 2454
regulate no food 2.38 2736
regulate yes autonomy 2.28 3279
regulate yes control 2.25 3549
regulate yes food 2.40 3780
NA yes autonomy 3.02 228
NA yes control 2.60 183
NA yes food 2.76 240

individual differences

ggplot(plot.rating, aes(action, meanRating, color = choice)) +
  geom_point(aes(group =  interaction(subjectID, choice), color = choice)) +
  geom_line(aes(group = interaction(subjectID, choice), color = choice)) + 
  facet_wrap(~condition + subjectID, ncol = 5) +
  geom_text(aes(x = 2.25, y = 3.5, label = taskChoice), color = "black", size = 3) + 
  scale_color_manual(values = c(palette[1],palette[3])) +
  labs(y = "mean craving rating") + 
  theme_minimal()
## Warning: Removed 268 rows containing missing values (geom_text).

difference in pre-task and task craving ratings by task condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition) %>%
  summarize(meanRating = mean(ratingDiff, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanRating)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = c(palette_cond, palette[2],palette[4])) +
  labs(y = "difference in craving rating") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanRating)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = c(palette_cond, palette[1],palette[3])) +
  labs(y = "difference in craving rating") + 
  theme_minimal()

data1 %>%
  group_by(action, choice) %>%
  summarize(meanRating = mean(ratingDiff, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean difference in craving ratings by condition")
mean difference in craving ratings by condition
action choice meanRating n
look no -0.52 7452
look yes -0.50 11061
regulate no -1.14 7428
regulate yes -1.15 10608
NA yes -0.66 651

difference in pre-task and task craving ratings by task and experimental condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition) %>%
  summarize(meanRating = mean(ratingDiff, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = palette_cond) +
  labs(y = "difference in craving rating") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = palette_cond) +
  labs(y = "difference in craving rating") + 
  theme_minimal()

data1 %>%
  group_by(action, choice, condition) %>%
  summarize(meanRating = mean(ratingDiff, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean difference in craving ratings by condition")
mean difference in craving ratings by condition
action choice condition meanRating n
look no autonomy -0.43 2262
look no control -0.56 2454
look no food -0.56 2736
look yes autonomy -0.41 3243
look yes control -0.52 3630
look yes food -0.54 4188
regulate no autonomy -1.05 2238
regulate no control -1.28 2454
regulate no food -1.09 2736
regulate yes autonomy -1.16 3279
regulate yes control -1.25 3549
regulate yes food -1.04 3780
NA yes autonomy -0.53 228
NA yes control -0.87 183
NA yes food -0.62 240

pre-task and task craving ratings by task and experimental condition

plot.rating = data1 %>%
  filter(!is.na(action)) %>%
  gather(ratingCondition, rating, respRating, likingRating) %>%
  mutate(ratingCondition = ifelse(ratingCondition == "respRating", "task rating",
                     ifelse(ratingCondition == "likingRating", "pre rating", condition))) %>%
  group_by(subjectID, action, choice, ratingCondition, condition) %>%
  summarize(meanRating = mean(rating, na.rm = TRUE))

ggplot(plot.rating, aes(action, meanRating)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(choice ~ ratingCondition) +
  scale_color_manual(values = c(palette_cond, palette[2],palette[4])) +
  labs(y = "mean craving rating") +
  theme_minimal()

mean craving ratings across trials

# linear trend
ggplot(filter(data1, !is.na(action)), aes(trial, respRating, color = action, linetype = choice)) + 
  geom_smooth(method = 'lm', se=FALSE, aes(group = interaction(subjectID, choice, action)), alpha =.06, size=.1) +
  geom_smooth(method = 'lm', se=FALSE, size = 1.5) + 
  facet_grid(~condition) + 
  scale_color_manual(values = c(palette[2],palette[4])) +
  labs(y = "mean craving rating") +
  theme_minimal()
## Warning: Removed 300 rows containing non-finite values (stat_smooth).

## Warning: Removed 300 rows containing non-finite values (stat_smooth).

mean craving rating by run

Visualizing the decrease over time by collapsing across run

plot.rating = data1 %>%
  group_by(subjectID, action, choice, run, condition) %>%
  summarize(meanRating = mean(respRating, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(run, meanRating, color = action, linetype = choice)) +
  geom_point(aes(group = interaction(subjectID, action, choice)), alpha = .05, size = 3) + 
  geom_line(aes(group = interaction(subjectID, action, choice)), alpha = .05, size = 1) + 
  stat_summary(aes(group = interaction(action, choice)), fun.data = "mean_cl_boot", size = 1) +
  stat_summary(aes(group = interaction(action, choice)), fun.y = "mean", geom = "line", size = 1) +
  facet_grid(~ condition) +
  scale_color_manual(values = c(palette[2], palette[4])) +
  labs(y = "mean craving rating") +
  theme_minimal()

data1 %>%
  group_by(action, choice, run, condition) %>%
  summarize(meanRating = mean(respRating, na.rm=TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean craving ratings by condition and run")
mean craving ratings by condition and run
action choice run condition meanRating n
look no run1 autonomy 3.00 762
look no run1 control 3.10 822
look no run1 food 3.01 924
look no run2 autonomy 3.13 750
look no run2 control 2.92 816
look no run2 food 2.85 900
look no run3 autonomy 2.98 750
look no run3 control 2.88 816
look no run3 food 2.86 912
look yes run1 autonomy 3.09 1053
look yes run1 control 3.09 1254
look yes run1 food 2.99 1380
look yes run2 autonomy 3.03 1041
look yes run2 control 2.94 1188
look yes run2 food 2.97 1392
look yes run3 autonomy 3.03 1149
look yes run3 control 2.94 1188
look yes run3 food 2.89 1416
regulate no run1 autonomy 2.40 738
regulate no run1 control 2.36 822
regulate no run1 food 2.49 900
regulate no run2 autonomy 2.44 750
regulate no run2 control 2.19 816
regulate no run2 food 2.34 924
regulate no run3 autonomy 2.38 750
regulate no run3 control 2.23 816
regulate no run3 food 2.30 912
regulate yes run1 autonomy 2.45 1101
regulate yes run1 control 2.34 1125
regulate yes run1 food 2.53 1248
regulate yes run2 autonomy 2.26 1125
regulate yes run2 control 2.19 1236
regulate yes run2 food 2.39 1296
regulate yes run3 autonomy 2.15 1053
regulate yes run3 control 2.22 1188
regulate yes run3 food 2.29 1236
NA yes run1 autonomy 2.75 96
NA yes run1 control 2.79 87
NA yes run1 food 2.92 108
NA yes run2 autonomy 3.19 84
NA yes run2 control 2.17 24
NA yes run2 food 2.00 48
NA yes run3 autonomy 3.30 48
NA yes run3 control 2.50 72
NA yes run3 food 3.00 84

test craving rating models

  • Model 1 = respRating ~ action*choice + (1 + action | subjectID)
  • Model 2 = respRating ~ action*choice*condition + (1 + action | subjectID)
  • Model 3 = respRating ~ action*choice*condition*taskChoice + (1 + action | subjectID)
# get complete data
data.complete = data1 %>%
  ungroup() %>%
  mutate(taskChoice = ifelse(is.na(taskChoice), "none", taskChoice),
         taskChoice = as.factor(taskChoice),
         condition = as.factor(condition)) %>%
  select(subjectID, condition, taskChoice, run, action, choice, trial, category, likingRating, craved, respRating, rtRating, ratingDiff) %>%
  na.omit()

# relevel condition
data.complete$condition = relevel(data.complete$condition, "control")

# run and compare models
model.1 = lmer(ratingDiff ~ action*choice + (1 + action | subjectID), 
               REML = FALSE,
               data = data.complete)
model.2 = lmer(ratingDiff ~ action*choice*condition + (1 + action | subjectID), 
               REML = FALSE,
               data = data.complete)
model.3 = lmer(ratingDiff ~ action*choice*taskChoice*condition + (1 + action | subjectID), 
               REML = FALSE,
               data = data.complete)

anova(model.1, model.2, model.3)
summary(model.1)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: ratingDiff ~ action * choice + (1 + action | subjectID)
##    Data: data.complete
## 
##      AIC      BIC   logLik deviance df.resid 
##  94662.7  94730.7 -47323.3  94646.7    36241 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9106 -0.6456  0.0678  0.6633  3.4209 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1648   0.4060        
##            actionregulate 0.2562   0.5062   -0.38
##  Residual                 0.7799   0.8831        
## Number of obs: 36249, groups:  subjectID, 105
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)              -5.217e-01  4.098e-02  1.138e+02 -12.730   <2e-16
## actionregulate           -6.180e-01  5.157e-02  1.155e+02 -11.984   <2e-16
## choiceyes                 2.024e-02  1.334e-02  3.606e+04   1.517    0.129
## actionregulate:choiceyes -2.207e-02  1.895e-02  3.605e+04  -1.165    0.244
##                             
## (Intercept)              ***
## actionregulate           ***
## choiceyes                   
## actionregulate:choiceyes    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys
## actionreglt -0.407              
## choiceyes   -0.193  0.153       
## actnrglt:ch  0.136 -0.216 -0.704
summary(model.2)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: 
## ratingDiff ~ action * choice * condition + (1 + action | subjectID)
##    Data: data.complete
## 
##      AIC      BIC   logLik deviance df.resid 
##    94644    94780   -47306    94612    36233 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9130 -0.6487  0.0710  0.6580  3.4170 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1618   0.4022        
##            actionregulate 0.2460   0.4960   -0.38
##  Residual                 0.7793   0.8828        
## Number of obs: 36249, groups:  subjectID, 105
## 
## Fixed effects:
##                                              Estimate Std. Error
## (Intercept)                                -5.588e-01  7.052e-02
## actionregulate                             -7.049e-01  8.787e-02
## choiceyes                                   3.020e-02  2.322e-02
## conditionautonomy                           1.222e-01  1.019e-01
## conditionfood                              -3.623e-04  9.756e-02
## actionregulate:choiceyes                   -2.114e-02  3.287e-02
## actionregulate:conditionautonomy            8.033e-02  1.270e-01
## actionregulate:conditionfood                1.732e-01  1.215e-01
## choiceyes:conditionautonomy                -1.364e-02  3.371e-02
## choiceyes:conditionfood                    -1.567e-02  3.189e-02
## actionregulate:choiceyes:conditionautonomy -9.823e-02  4.768e-02
## actionregulate:choiceyes:conditionfood      7.983e-02  4.538e-02
##                                                    df t value Pr(>|t|)    
## (Intercept)                                 1.147e+02  -7.924 1.65e-12 ***
## actionregulate                              1.165e+02  -8.023 9.09e-13 ***
## choiceyes                                   3.605e+04   1.301   0.1934    
## conditionautonomy                           1.143e+02   1.199   0.2330    
## conditionfood                               1.140e+02  -0.004   0.9970    
## actionregulate:choiceyes                    3.605e+04  -0.643   0.5201    
## actionregulate:conditionautonomy            1.162e+02   0.633   0.5282    
## actionregulate:conditionfood                1.158e+02   1.425   0.1569    
## choiceyes:conditionautonomy                 3.606e+04  -0.405   0.6857    
## choiceyes:conditionfood                     3.605e+04  -0.491   0.6232    
## actionregulate:choiceyes:conditionautonomy  3.605e+04  -2.060   0.0394 *  
## actionregulate:choiceyes:conditionfood      3.605e+04   1.759   0.0786 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                       (Intr) actnrg chocys cndtnt cndtnf actnrglt:ch
## actionreglt           -0.407                                        
## choiceyes             -0.195  0.157                                 
## conditntnmy           -0.692  0.282  0.135                          
## conditionfd           -0.723  0.294  0.141  0.500                   
## actnrglt:ch            0.138 -0.221 -0.706 -0.095 -0.100            
## actnrglt:cndtnt        0.282 -0.692 -0.108 -0.406 -0.204  0.153     
## actnrglt:cndtnf        0.294 -0.723 -0.113 -0.204 -0.406  0.160     
## chcys:cndtnt           0.134 -0.108 -0.689 -0.194 -0.097  0.487     
## chcys:cndtnf           0.142 -0.114 -0.728 -0.098 -0.195  0.514     
## actnrglt:chcys:cndtnt -0.095  0.152  0.487  0.137  0.069 -0.689     
## actnrglt:chcys:cndtnf -0.100  0.160  0.512  0.069  0.137 -0.724     
##                       actnrglt:cndtnt actnrglt:cndtnf chcys:cndtnt
## actionreglt                                                       
## choiceyes                                                         
## conditntnmy                                                       
## conditionfd                                                       
## actnrglt:ch                                                       
## actnrglt:cndtnt                                                   
## actnrglt:cndtnf        0.500                                      
## chcys:cndtnt           0.156           0.078                      
## chcys:cndtnf           0.079           0.157           0.502      
## actnrglt:chcys:cndtnt -0.221          -0.110          -0.707      
## actnrglt:chcys:cndtnf -0.111          -0.220          -0.352      
##                       chcys:cndtnf actnrglt:chcys:cndtnt
## actionreglt                                             
## choiceyes                                               
## conditntnmy                                             
## conditionfd                                             
## actnrglt:ch                                             
## actnrglt:cndtnt                                         
## actnrglt:cndtnf                                         
## chcys:cndtnt                                            
## chcys:cndtnf                                            
## actnrglt:chcys:cndtnt -0.355                            
## actnrglt:chcys:cndtnf -0.703        0.499
# summarize
control = lmer(ratingDiff ~ action*choice + (1 + action | subjectID), 
               REML = FALSE,
               data = filter(data.complete, condition == "control"))
autonomy = lmer(ratingDiff ~ action*choice + (1 + action | subjectID), 
                REML = FALSE, 
               data = filter(data.complete, condition == "autonomy"))
food = lmer(ratingDiff ~ action*choice*taskChoice + (1 + action | subjectID), 
            REML = FALSE,
            data = filter(data.complete, condition == "food"))

summary(control)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: ratingDiff ~ action * choice + (1 + action | subjectID)
##    Data: filter(data.complete, condition == "control")
## 
##      AIC      BIC   logLik deviance df.resid 
##  29552.5  29611.7 -14768.3  29536.5    12015 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6334 -0.7191  0.0987  0.6163  3.3805 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1494   0.3865        
##            actionregulate 0.2375   0.4873   -0.45
##  Residual                 0.6683   0.8175        
## Number of obs: 12023, groups:  subjectID, 35
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)              -5.590e-01  6.760e-02  3.799e+01  -8.270 5.08e-10
## actionregulate           -7.046e-01  8.591e-02  3.816e+01  -8.202 6.02e-10
## choiceyes                 3.015e-02  2.150e-02  1.196e+04   1.402    0.161
## actionregulate:choiceyes -2.112e-02  3.044e-02  1.196e+04  -0.694    0.488
##                             
## (Intercept)              ***
## actionregulate           ***
## choiceyes                   
## actionregulate:choiceyes    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys
## actionreglt -0.470              
## choiceyes   -0.188  0.148       
## actnrglt:ch  0.133 -0.209 -0.706
summary(autonomy)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: ratingDiff ~ action * choice + (1 + action | subjectID)
##    Data: filter(data.complete, condition == "autonomy")
## 
##      AIC      BIC   logLik deviance df.resid 
##  30145.0  30203.4 -15064.5  30129.0    10894 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6241 -0.6414  0.0920  0.6023  3.1646 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1786   0.4226        
##            actionregulate 0.2598   0.5097   -0.27
##  Residual                 0.9087   0.9532        
## Number of obs: 10902, groups:  subjectID, 32
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)              -4.367e-01  7.744e-02  3.495e+01  -5.640 2.32e-06
## actionregulate           -6.244e-01  9.462e-02  3.559e+01  -6.599 1.17e-07
## choiceyes                 1.660e-02  2.639e-02  1.085e+04   0.629  0.52924
## actionregulate:choiceyes -1.194e-01  3.729e-02  1.084e+04  -3.202  0.00137
##                             
## (Intercept)              ***
## actionregulate           ***
## choiceyes                   
## actionregulate:choiceyes ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys
## actionreglt -0.306              
## choiceyes   -0.199  0.163       
## actnrglt:ch  0.141 -0.231 -0.708
summary(food)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: 
## ratingDiff ~ action * choice * taskChoice + (1 + action | subjectID)
##    Data: filter(data.complete, condition == "food")
## 
##      AIC      BIC   logLik deviance df.resid 
##  34697.8  34787.8 -17336.9  34673.8    13312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7115 -0.5712 -0.0017  0.7085  3.2550 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1588   0.3985        
##            actionregulate 0.2409   0.4908   -0.42
##  Residual                 0.7736   0.8795        
## Number of obs: 13324, groups:  subjectID, 38
## 
## Fixed effects:
##                                               Estimate Std. Error
## (Intercept)                                 -5.633e-01  8.781e-02
## actionregulate                              -5.196e-01  1.093e-01
## choiceyes                                    8.341e-03  2.876e-02
## taskChoiceregulate                           9.813e-03  1.353e-01
## actionregulate:choiceyes                     8.894e-02  4.099e-02
## actionregulate:taskChoiceregulate           -2.905e-02  1.684e-01
## choiceyes:taskChoiceregulate                 1.446e-02  4.402e-02
## actionregulate:choiceyes:taskChoiceregulate -7.240e-02  6.314e-02
##                                                     df t value Pr(>|t|)
## (Intercept)                                  4.102e+01  -6.414 1.11e-07
## actionregulate                               4.204e+01  -4.755 2.34e-05
## choiceyes                                    1.325e+04   0.290    0.772
## taskChoiceregulate                           4.101e+01   0.073    0.943
## actionregulate:choiceyes                     1.325e+04   2.170    0.030
## actionregulate:taskChoiceregulate            4.202e+01  -0.173    0.864
## choiceyes:taskChoiceregulate                 1.325e+04   0.329    0.742
## actionregulate:choiceyes:taskChoiceregulate  1.325e+04  -1.147    0.252
##                                                
## (Intercept)                                 ***
## actionregulate                              ***
## choiceyes                                      
## taskChoiceregulate                             
## actionregulate:choiceyes                    *  
## actionregulate:taskChoiceregulate              
## choiceyes:taskChoiceregulate                   
## actionregulate:choiceyes:taskChoiceregulate    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys tskChc actnr: actn:C chcy:C
## actionreglt -0.443                                          
## choiceyes   -0.195  0.157                                   
## taskChcrglt -0.649  0.287  0.127                            
## actnrglt:ch  0.137 -0.221 -0.702 -0.089                     
## actnrglt:tC  0.287 -0.649 -0.102 -0.443  0.143              
## chcys:tskCh  0.128 -0.103 -0.653 -0.196  0.458  0.158       
## actnrglt::C -0.089  0.143  0.456  0.137 -0.649 -0.221 -0.697

effort ratings

by subject and task condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition) %>%
  summarize(meanEffort = mean(respEffort, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanEffort)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = c(palette_cond, palette[2],palette[4])) +
  labs(y = "mean effort rating") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanEffort)) +
  geom_point(aes(group = subjectID, color = condition), alpha = .05, size = 3) + 
  geom_line(aes(group = subjectID, color = condition), alpha = .05, size = 1) + 
  stat_summary(aes(group = 1), color = palette_grey, fun.y = mean, geom = "line") +
  stat_summary(color = palette_grey, fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = c(palette_cond, palette[1],palette[3])) +
  labs(y = "mean effort rating") + 
  theme_minimal()

data1 %>%
  group_by(action, choice) %>%
  summarize(meanEffort = mean(respEffort, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean effort ratings by condition")
mean effort ratings by condition
action choice meanEffort n
look no 1.68 7452
look yes 1.66 11061
regulate no 2.10 7428
regulate yes 2.09 10608
NA yes 1.81 651

by task and experimental condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition, taskChoice) %>%
  summarize(meanEffort = mean(respEffort, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanEffort, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean effort rating") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanEffort, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean effort rating") + 
  theme_minimal()

data1 %>%
  group_by(action, choice, condition) %>%
  summarize(meanEffort = mean(respEffort, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean effort ratings by condition")
mean effort ratings by condition
action choice condition meanEffort n
look no autonomy 1.61 2262
look no control 1.74 2454
look no food 1.70 2736
look yes autonomy 1.60 3243
look yes control 1.72 3630
look yes food 1.65 4188
regulate no autonomy 1.99 2238
regulate no control 2.15 2454
regulate no food 2.15 2736
regulate yes autonomy 1.96 3279
regulate yes control 2.13 3549
regulate yes food 2.17 3780
NA yes autonomy 1.73 228
NA yes control 2.07 183
NA yes food 1.68 240

individual differences

ggplot(plot.rating, aes(action, meanEffort, color = choice)) +
  geom_point(aes(group =  interaction(subjectID, choice), color = choice)) +
  geom_line(aes(group = interaction(subjectID, choice), color = choice)) + 
  facet_wrap(~condition + subjectID, ncol = 5) +
  geom_text(aes(x = 2.25, y = 3.5, label = taskChoice), color = "black", size = 3) + 
  scale_color_manual(values = c(palette[1],palette[3])) +
  labs(y = "mean effort rating") + 
  theme_minimal()
## Warning: Removed 268 rows containing missing values (geom_text).

test effort rating models

  • Model 1 = respEffort ~ action*choice + (1 + action | subjectID)
  • Model 2 = respEffort ~ action*choice*condition + (1 + action | subjectID)
  • Model 3 = respEffort ~ action*choice*condition*taskChoice + (1 + action | subjectID)
# get complete data
data.complete = data1 %>%
  ungroup() %>%
  mutate(taskChoice = ifelse(is.na(taskChoice), "none", taskChoice),
         taskChoice = as.factor(taskChoice),
         condition = as.factor(condition)) %>%
  select(subjectID, condition, taskChoice, run, action, choice, trial, category, likingRating, craved, respEffort, rtEffort) %>%
  na.omit()

# relevel condition
data.complete$condition = relevel(data.complete$condition, "control")

# run and compare models
model.1 = lmer(respEffort ~ action*choice + (1 + action | subjectID), 
               REML = FALSE,
               data = data.complete)
model.2 = lmer(respEffort ~ action*choice*condition + (1 + action | subjectID), 
               REML = FALSE,
               data = data.complete)
model.3 = lmer(respEffort ~ action*choice*condition*taskChoice + (1 + action | subjectID), 
               REML = FALSE,
               data = data.complete)

anova(model.1, model.2, model.3)
summary(model.1)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: respEffort ~ action * choice + (1 + action | subjectID)
##    Data: data.complete
## 
##      AIC      BIC   logLik deviance df.resid 
##  89773.3  89841.1 -44878.7  89757.3    35497 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8133 -0.7316 -0.1783  0.6728  3.4265 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1265   0.3557        
##            actionregulate 0.2210   0.4701   -0.40
##  Residual                 0.7181   0.8474        
## Number of obs: 35505, groups:  subjectID, 105
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)               1.689e+00  3.617e-02  1.149e+02  46.697  < 2e-16
## actionregulate            4.062e-01  4.806e-02  1.164e+02   8.453 9.41e-14
## choiceyes                -1.721e-02  1.294e-02  3.532e+04  -1.330    0.183
## actionregulate:choiceyes -2.780e-03  1.838e-02  3.531e+04  -0.151    0.880
##                             
## (Intercept)              ***
## actionregulate           ***
## choiceyes                   
## actionregulate:choiceyes    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys
## actionreglt -0.424              
## choiceyes   -0.212  0.159       
## actnrglt:ch  0.149 -0.225 -0.704
# summarize
control = lmer(respEffort ~ action*choice + (1 + action | subjectID), 
               REML = FALSE,
               data = filter(data.complete, condition == "control"))
autonomy = lmer(respEffort ~ action*choice + (1 + action | subjectID), 
                REML = FALSE, 
               data = filter(data.complete, condition == "autonomy"))
food = lmer(respEffort ~ action*choice*taskChoice + (1 + action | subjectID), 
            REML = FALSE,
            data = filter(data.complete, condition == "food"))

summary(control)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: respEffort ~ action * choice + (1 + action | subjectID)
##    Data: filter(data.complete, condition == "control")
## 
##      AIC      BIC   logLik deviance df.resid 
##  29545.8  29604.8 -14764.9  29529.8    11783 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1838 -0.6967 -0.1835  0.7142  3.1385 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1622   0.4027        
##            actionregulate 0.2201   0.4691   -0.64
##  Residual                 0.7016   0.8376        
## Number of obs: 11791, groups:  subjectID, 35
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)               1.748e+00  7.040e-02  3.755e+01  24.835  < 2e-16
## actionregulate            4.014e-01  8.320e-02  3.915e+01   4.824 2.16e-05
## choiceyes                 2.627e-03  2.223e-02  1.173e+04   0.118    0.906
## actionregulate:choiceyes -2.114e-02  3.150e-02  1.172e+04  -0.671    0.502
##                             
## (Intercept)              ***
## actionregulate           ***
## choiceyes                   
## actionregulate:choiceyes    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys
## actionreglt -0.647              
## choiceyes   -0.186  0.158       
## actnrglt:ch  0.132 -0.224 -0.706
summary(autonomy)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: respEffort ~ action * choice + (1 + action | subjectID)
##    Data: filter(data.complete, condition == "autonomy")
## 
##      AIC      BIC   logLik deviance df.resid 
##  26173.4  26231.6 -13078.7  26157.4    10682 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9293 -0.7150 -0.2337  0.6287  3.3491 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.1210   0.3479        
##            actionregulate 0.2376   0.4874   -0.26
##  Residual                 0.6615   0.8133        
## Number of obs: 10690, groups:  subjectID, 32
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)               1.609e+00  6.395e-02  3.505e+01  25.158  < 2e-16
## actionregulate            3.611e-01  8.967e-02  3.476e+01   4.027 0.000291
## choiceyes                -5.124e-03  2.274e-02  1.063e+04  -0.225 0.821732
## actionregulate:choiceyes -3.160e-02  3.211e-02  1.063e+04  -0.984 0.325212
##                             
## (Intercept)              ***
## actionregulate           ***
## choiceyes                   
## actionregulate:choiceyes    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys
## actionreglt -0.289              
## choiceyes   -0.207  0.148       
## actnrglt:ch  0.146 -0.209 -0.708
summary(food)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: 
## respEffort ~ action * choice * taskChoice + (1 + action | subjectID)
##    Data: filter(data.complete, condition == "food")
## 
##      AIC      BIC   logLik deviance df.resid 
##  33975.6  34065.3 -16975.8  33951.6    13012 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3164 -0.7665 -0.1445  0.6717  3.2427 
## 
## Random effects:
##  Groups    Name           Variance Std.Dev. Corr 
##  subjectID (Intercept)    0.08416  0.2901        
##            actionregulate 0.20149  0.4489   -0.30
##  Residual                 0.77839  0.8823        
## Number of obs: 13024, groups:  subjectID, 38
## 
## Fixed effects:
##                                               Estimate Std. Error
## (Intercept)                                  1.654e+00  6.582e-02
## actionregulate                               3.919e-01  1.008e-01
## choiceyes                                   -7.679e-02  2.916e-02
## taskChoiceregulate                           1.144e-01  1.014e-01
## actionregulate:choiceyes                     1.409e-01  4.149e-02
## actionregulate:taskChoiceregulate            1.315e-01  1.554e-01
## choiceyes:taskChoiceregulate                 7.581e-02  4.474e-02
## actionregulate:choiceyes:taskChoiceregulate -2.478e-01  6.413e-02
##                                                     df t value Pr(>|t|)
## (Intercept)                                  4.406e+01  25.126  < 2e-16
## actionregulate                               4.286e+01   3.886 0.000349
## choiceyes                                    1.296e+04  -2.634 0.008457
## taskChoiceregulate                           4.409e+01   1.128 0.265430
## actionregulate:choiceyes                     1.296e+04   3.394 0.000690
## actionregulate:taskChoiceregulate            4.287e+01   0.846 0.402279
## choiceyes:taskChoiceregulate                 1.296e+04   1.695 0.090194
## actionregulate:choiceyes:taskChoiceregulate  1.296e+04  -3.865 0.000112
##                                                
## (Intercept)                                 ***
## actionregulate                              ***
## choiceyes                                   ** 
## taskChoiceregulate                             
## actionregulate:choiceyes                    ***
## actionregulate:taskChoiceregulate              
## choiceyes:taskChoiceregulate                .  
## actionregulate:choiceyes:taskChoiceregulate ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) actnrg chocys tskChc actnr: actn:C chcy:C
## actionreglt -0.344                                          
## choiceyes   -0.264  0.172                                   
## taskChcrglt -0.649  0.223  0.171                            
## actnrglt:ch  0.186 -0.242 -0.703 -0.120                     
## actnrglt:tC  0.223 -0.649 -0.112 -0.344  0.157              
## chcys:tskCh  0.172 -0.112 -0.652 -0.266  0.458  0.174       
## actnrglt::C -0.120  0.157  0.455  0.186 -0.647 -0.242 -0.698

reaction times

craving RTs by condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition, taskChoice) %>%
  summarize(meanRating = mean(rtRating, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean craving RT") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean craving RT") + 
  theme_minimal()

data1 %>%
  group_by(action, choice, condition) %>%
  summarize(meanEffort = mean(respRating, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean effort RTs by condition")
mean effort RTs by condition
action choice condition meanEffort n
look no autonomy 3.04 2262
look no control 2.97 2454
look no food 2.91 2736
look yes autonomy 3.05 3243
look yes control 2.99 3630
look yes food 2.95 4188
regulate no autonomy 2.40 2238
regulate no control 2.26 2454
regulate no food 2.38 2736
regulate yes autonomy 2.28 3279
regulate yes control 2.25 3549
regulate yes food 2.40 3780
NA yes autonomy 3.02 228
NA yes control 2.60 183
NA yes food 2.76 240

craving RTs by condition and rating

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition, respRating, taskChoice) %>%
  summarize(meanRating = mean(rtRating, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(filter(plot.rating, !is.na(respRating)), aes(action, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(choice~respRating) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean craving RT") + 
  theme_minimal()

ggplot(filter(plot.rating, !is.na(respRating)), aes(choice, meanRating, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(action~respRating) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean craving RT") + 
  theme_minimal()

data1 %>%
  group_by(action, choice, condition) %>%
  summarize(meanEffort = mean(respRating, na.rm = TRUE),
            n = n()) %>%
  kable(digits = 2, format = "pandoc", caption = "mean effort RTs by condition")
mean effort RTs by condition
action choice condition meanEffort n
look no autonomy 3.04 2262
look no control 2.97 2454
look no food 2.91 2736
look yes autonomy 3.05 3243
look yes control 2.99 3630
look yes food 2.95 4188
regulate no autonomy 2.40 2238
regulate no control 2.26 2454
regulate no food 2.38 2736
regulate yes autonomy 2.28 3279
regulate yes control 2.25 3549
regulate yes food 2.40 3780
NA yes autonomy 3.02 228
NA yes control 2.60 183
NA yes food 2.76 240

effort RTs by condition

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition, taskChoice) %>%
  summarize(meanEffort = mean(rtEffort, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(plot.rating, aes(action, meanEffort, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~choice) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean effort RT") + 
  theme_minimal()

ggplot(plot.rating, aes(choice, meanEffort, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(~action) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean effort RT") + 
  theme_minimal()

effort RTs by condition and rating

plot.rating = data1 %>%
  group_by(subjectID, action, choice, condition, respEffort, taskChoice) %>%
  summarize(meanEffort = mean(rtEffort, na.rm = TRUE)) %>%
  filter(!is.na(action))

ggplot(filter(plot.rating, !is.na(respEffort)), aes(action, meanEffort, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(choice~respEffort) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean effort RT") + 
  theme_minimal()

ggplot(filter(plot.rating, !is.na(respEffort)), aes(choice, meanEffort, color = condition)) +
  stat_summary(aes(group = condition), fun.y = mean, geom = "line") +
  stat_summary(aes(color = condition), fun.data = "mean_cl_boot", size = 1) +
  facet_grid(action~respEffort) +
  scale_color_manual(values = palette_cond) +
  labs(y = "mean effort RT") + 
  theme_minimal()

liking ratings (pre-task)

ratings by category and frequency chosen

These are only the ratings for the 90 images participants saw during the tasks

plot.rating = data1 %>%
  group_by(category) %>%
  summarize(meanRating = mean(likingRating, na.rm = TRUE),
            n = n(),
            sd = sd(likingRating, na.rm = TRUE),
            se = sd/sqrt(n),
            upper = meanRating + 2*se,
            lower = meanRating - 2*se)

ggplot(plot.rating, aes(reorder(category, -n), meanRating, color = category)) + 
  geom_point(aes(size = n)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + 
  scale_color_manual(values = wesanderson::wes_palette("Zissou1", 15, type = "continuous")) +
  labs(x = "category", y = "mean liking rating") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

distribution by subject and category

ggplot(data1, aes(likingRating, fill = category)) + 
  geom_histogram(stat = "bin", binwidth = .5) + 
  facet_wrap(~subjectID, ncol = 5) +
  scale_fill_manual(values = wesanderson::wes_palette("Zissou1", 15, type = "continuous")) +
  labs(x = "liking rating") +
  theme_minimal()

identify problematic categories and images

These are the ratings for all 45 images in each category

load rating task data

# define variables and paths
sub_dir = "~/Dropbox (PfeiBer Lab)/FreshmanProject/Tasks/ROC-C/output/"
sub_pattern = "HL[0-9]{3}"
subjects = list.files(sub_dir, pattern = sub_pattern)
runs = c("ratings")

# initialize data frame
ratings = data.frame()

# load data
for (sub in subjects) {
  for (run in runs) {
    file = paste0(sub_dir, sub, '/', sub, '_', run,'.csv')
    tmp = tryCatch(read.csv(file, stringsAsFactors = FALSE, header = FALSE) %>%
                     rename("rating" = V1,
                            "ranking" = V2,
                            "image" = V3) %>%
                     mutate(subjectID = sub) %>%
                     extract(image, "category", "(.*)[0-9]{2}.jpg", remove = FALSE), error = function(e) NULL)
      ratings = bind_rows(ratings, tmp)
      rm(tmp)
  }
}

ratings by category and frequency chosen

plot.rating = ratings %>%
  group_by(category) %>%
  summarize(meanRating = mean(rating, na.rm = TRUE),
            n = n(),
            sd = sd(rating, na.rm = TRUE),
            se = sd/sqrt(n),
            upper = meanRating + 2*se,
            lower = meanRating - 2*se)

ggplot(plot.rating, aes(reorder(category, -n), meanRating, color = category)) + 
  geom_point(aes(size = n)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + 
  scale_color_manual(values = wesanderson::wes_palette("Zissou1", 15, type = "continuous")) +
  labs(x = "category", y = "mean liking rating") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

summarize data

mean.ratings = ratings %>%
  group_by(image, category) %>%
  summarize(meanRating = mean(rating, na.rm = T),
            nChosen = n()) %>%
  arrange(desc(meanRating))

# higher than 2.5
mean.ratings %>%
  filter(meanRating > 2.5)
# lower than 2.5
mean.ratings %>%
  filter(meanRating <= 2.5) 
# good categories
mean.ratings %>%
  filter(meanRating > 2.5) %>%
  group_by(category, nChosen) %>%
  summarize(n = n()) %>%
  arrange(desc(n))
# problematic categories
mean.ratings %>%
  filter(meanRating <= 2.5) %>%
  group_by(category, nChosen) %>%
  summarize(n = n()) %>%
  arrange(desc(n))